#include "Temp.h"
#include "Global.h"
#include <algorithm>
#include <iomanip>
#include <iostream>
#include <cmath>
#include "FilesManager.h"
#include "TS.h"
#include <algorithm>

#ifdef GPU_INCLUDED
        #include "GPGPU.h"
#endif

Temp::Temp()
{
    static_.wyzn = 0.0;
    epsilon = std::numeric_limits<double>::epsilon();
    temhist.iidt = 1;
    temhist.ildt = 1;
}


//cl::Buffer *A = new cl::Buffer(*GPU.context, CL_MEM_READ_WRITE, sizeof (double) * ((NWMAX * MWMAX)*(NWMAX + 2)));

#ifdef GPU_INCLUDED
    cl::Event event;
#endif

int Temp::gauss1(int nn, int mm, double *a)
{

#ifdef GPU_INCLUDED   

         try
    {
    (*GPU.queue).enqueueWriteBuffer(*GPU.A, CL_TRUE, 0, sizeof (double) * ((NWMAX * MWMAX)*(NWMAX + 2)), a);//, NULL, &event);

    (*GPU.kernel_gauss).setArg(0, *GPU.A);
    (*GPU.kernel_gauss).setArg(1, sizeof (int), &nn);
    (*GPU.kernel_gauss).setArg(2, sizeof (int), &mm);

    for (int n = 0; n < nn - 1; n++)
    {

        (*GPU.kernel_gauss).setArg(3, sizeof (int), &n);

        (*GPU.queue).enqueueNDRangeKernel((*GPU.kernel_gauss), cl::NullRange, cl::NDRange(512), cl::NullRange);//, NULL, &event);

    }

    (*GPU.queue).enqueueReadBuffer(*GPU.A, CL_TRUE, 0, sizeof (double) * ((NWMAX * MWMAX)*(NWMAX + 2)), a);//, NULL, &event);
}
    catch (cl::Error& err)
    {
        std::cerr << "caught exception: " << err.what() << '(' << err.err() << ')' << std::endl;
        exit(1);
    }
#else
        
        double c;
       int k, l, mb;
        for (int n = 0; n < nn - 1; n++)
        {
            if (fabs(a[n]) > epsilon)
            {
                mb = std::min(mm, nn - n);
//                std::cout << "?MI " << mb << " " << std::flush;

                for (l = 0; l < mb-1; ++l)
                {
                    c = a[n + MWMAX * NWMAX * (l+1)] / a[n];

                    if (fabs(c) > epsilon)
                    {
                        for (k = 0; k < mb-l-1; ++k)
                        {
                            a[n + MWMAX * NWMAX * k + (l+1)] -= c * a[n + MWMAX * NWMAX * (k+l+1)];
                        }
                        
                        a[n + (l+1) * MWMAX * NWMAX] = c;   

                    }
                }
            }
        }
#endif


    return 0;

}

int Temp::gauss2(int nn, int mm, double *a, double *b)
{
    double c;
    int m, n, ih;

    for (n = 0; n < nn; n++)
    {
        c = b[n];
        if (std::abs(a[n]) > epsilon)
        {
            b[n] /= a[n];
        }

        ih = std::min(nn, n + mm);
        m = n + 1;
        for (int i = n + 1; i < ih; ++i)
        {
            m += MWMAX*NWMAX;
            b[i] -= a[m - 1] * c;
        }
    }

    for (int q = n - 1; q > 0; q--)
    {
        ih = std::min(nn, q + mm - 1);
        m = q;
        for (int i = q; i < ih; ++i)
        {
            m += MWMAX*NWMAX;
            b[q - 1] -= a[m - 1] * b[i];
        }
    }

    return 0;
}

void pokaz(double *qq, int iw, int npas)
{
    //(NWMAX + 2)*(MWMAX * NWMAX)
    //double *q = qq;
    //q -= MWMAX*NWMAX+1;
    std::cout << std::endl;
    std::cout << std::endl;
    std::cout << std::endl;
    std::cout << std::endl;
    /*  for (int i = 0; i < nw / (npas - 2); i++)
      {
          for (int j = 0; j < nw; j++)
              std::cout << qq[i * nw / (npas - 2) + j] << " ";
          std::cout << std::endl;
      }
     */
    for (int i = 0; i < iw; ++i)
    {
        //f[i] = 0.;
        for (int j = 0; j < npas; ++j)
        {
            //h[i + j * MWMAX * NWMAX] = 0.;
            std::cout << std::setprecision(2) << qq[i + j * MWMAX * NWMAX] << "\t";
            //  c[i + j * MWMAX * NWMAX] = 0.;
        }
        std::cout << std::endl;
    }
}

int Temp::temper(bool iosym, int nrchl, int mw, int nw, int ndt, double dt, int
                 iprint, int *idt, int istem, double temgrd, double *dtemgrd, int iwn, double tempgrd, double *dtempgrd)
{

    

    double dysyp = 1.; /*                  dysyp - wspolczynnik dysypacji energii odksztalcenia */
    int j1, j2;
    double cc;
    int ij, ji, it, mdt;
    int npas1;
    double tchcr;
    double xyzak;
    int ie = (mw - 1) * (nw - 1);
    int iw = mw * nw;
    //    int ip = iw << 1;
    int npas = nw + 2;
    double xmax = -1e10;
    double xmin = 1e10;

    for (int i = 0; i < iw; ++i)
    {
        xmax = std::max(GLOB.zargd.xw[2 * i], xmax);
        xmin = std::min(GLOB.zargd.xw[2 * i], xmin);
    }
    xyzak = xmax - xmin;


    AddTchcal(GLOB.ral.tchcal);
    //AddNdt(ndt);
    AddT(GLOB.tp);

    tchcr = 0.;
    mdt = 1;
    /* -------------------------------- */
    /*     mdt=10 */
    //    double mcont = 0.;
    if (nrchl > 0)
    {
        *idt = 0;
        mdt = 1;
        /* -------------------------------- */
        /*        mdt=10 */
        //      mcont = 0.;
        savetem(&mdt, idt, tchcr, mw, nw);
    }

    int ipr = 0;
    for (int ik = 0; ik < ndt; ++ik)
    {
        //std::cout << "przej " << ik << " z " << ndt << std::endl;
        /*        ik=ik+1 */
        it = ik + 1;

        if (ik + 1 == ndt)
        {
            ipr = iprint - 1;
        }

        for (int i = 0; i < iw; ++i)
        {
            f[i] = 0.;
            for (int j = 0; j < npas; ++j)
            {
                h[i + j * MWMAX * NWMAX] = 0.;
                c[i + j * MWMAX * NWMAX] = 0.;
            }
        }
        npas1 = 0;
        ij = 0;
        ji = 0;

        for (int k = 0; k < ie; ++k)
        {
            elem(iosym, k + 1, it, &iw, xyzak / 1000., &dysyp, dt, &npas1, &ij, &ji, nrchl);
        }
        if (npas1 > npas)
        {
            std::cout << "blad macierzy h  -  zalozona szerokosc " << npas1 << " " << npas << " zostala przekroczona" << std::endl;
            exit(1);
        }

        for (int i = 1; i <= iw; ++i)
        {
            cc = 0.;
            if (i != 1)
            {
                j1 = std::max(1, i - npas + 1);

                for (int j = j1; j <= i - 1; ++j)
                {
                    cc += c[j + (i - j + 1) * MWMAX * NWMAX - (MWMAX * NWMAX + 1)] * GLOB.tp[j - 1];
                }
            }

            j2 = std::min(npas + i - 1, iw);
            for (int j = i; j <= j2; ++j)
            {
                cc += c[i + (j - i + 1) * MWMAX * NWMAX - (MWMAX * NWMAX + 1)] * GLOB.tp[j - 1];
            }
            f[i - 1] = cc - f[i - 1] * 3.;
        }

        //     int it = 0;
        //   for (int q = 0; q < (NWMAX + 2); q++)
        //     for (int v = 0; v < (NWMAX * MWMAX); v++)
        //       hh[q + v * M] = h[it++];

        //double *nnq = h;
        //nnq -= MWMAX * NWMAX + 1;
        // pokaz(h, iw, npas);
        //std::cout << iw << " " << npas << " " << (NWMAX + 2)*(MWMAX * NWMAX) << std::endl;
        //gauss1(iw, npas, nnq);
        GLOB.LICZ++;
        gauss1(iw, npas, h);
//        cl_ulong start = event.getProfilingInfo<CL_PROFILING_COMMAND_START>();
//        cl_ulong end = event.getProfilingInfo<CL_PROFILING_COMMAND_END>();
//        std::cout << "Time " << (end-start) / 1000000.<< " ms " << std::endl;
//        gauss1(iw, npas, h);

        //it = 0;
        //for (int q = 0; q < (NWMAX + 2); q++)
        //  for (int v = 0; v < (NWMAX * MWMAX); v++)
        //    h[it++] = hh[q + v * M];


        //h = nnq;
        // pokaz(h, iw, npas);
        gauss2(iw, npas, h, f);


        // pokaz(h, iw, npas);
        //   exit(1);
        ++ipr;
        tchcr += dt;
        GLOB.ral.tchcal += dt;
        putem(istem, f, GLOB.tp, GLOB.ral.tchcal, GLOB.nstem, temgrd, dtemgrd);
        putemp(iwn, f, GLOB.ral.tchcal, iw, nw, tempgrd, dtempgrd);

        for (int i = 0; i < iw; ++i)
        {
            if (GLOB.tprzem[i] < 0.)
            {
                if (temhist.iidt > ILDTMAX)
                {
                    for (int j = 2; j <= temhist.ildt; ++j)
                    {
                        temhist.temhist[i][j - 2] = temhist.temhist[i][j - 1];
                    }
                }
                temhist.temhist[i][temhist.ildt - 1] = GLOB.tp[i];
                temhist.temlast[i] = f[i];
            }
            crate[i] = (temhist.temhist[i][0] - temhist.temlast[i]) / dt / temhist.ildt;
            if (crate[i] < 1e-10)
            {
                crate[i] = 1e-10;
            }
            /*            xv=xvv(i) */
            /*           Q(i)=cieprzem(dt,tp(i),Crate(i),tchcal,tprzem(i),xvv(i)) */
            if (f[i] >= GLOB.ral.tpow)
            {
                GLOB.tp[i] = f[i];
            }
            else
            {
                GLOB.tp[i] = GLOB.ral.tpow;
            }
            GLOB.q[i] = cieprzem(dt, GLOB.tp[i], crate[i], GLOB.ral.tchcal, &GLOB.tprzem[i], &GLOB.xvv[i]);
        }
        ++temhist.iidt;
        if (temhist.ildt < ILDTMAX)
        {
            ++temhist.ildt;
        }
        if (ipr >= iprint)
        {
            ipr = 0;
            if (ik == 0)
            {
                ipr = 1;
            }

            ShowTemps(true);
            AddTchcal(GLOB.ral.tchcal);
            AddT(GLOB.tp);
        }
        else
        {
            ShowTemps(false);
        }
    }

    AddTFinal(); // dodaj do v_tf tablice temperatur bedacych na szczycie tab. v_t

    if (nrchl > 0)
    {
        --(*idt);
    }

    return 0;

}

int Temp::savetem(int *mdt, int *idt, double tchcr, int mw, int nw)
{

    int k;

    if (*idt == 0)
    {
        GLOB.tpc[1][0] = tchcr;
        k = 1;
        for (int i = 1; i <= mw; ++i)
        {
            GLOB.tpc[0][i] = GLOB.zargd.xw[2 * k - 2] * 1000.;
            GLOB.tpc[1][i] = GLOB.tp[k - 1];
            k += nw;
        }
    }
    else
    {
        GLOB.tpc[*idt + 1][0] = tchcr;
        k = 1;
        for (int i = 1; i <= mw; ++i)
        {
            GLOB.tpc[*idt + 1][i] = GLOB.tp[k - 1];
            k += nw;
        }
    }
    ++(*idt);
    if (*idt > 10)
    {
        /* -------------------------------- */
        /*     if (idt.gt.25) then */
        for (int j = 1; j <= 5; ++j)
        {
            for (int i = 0; i < mw + 1; ++i)
            {
                GLOB.tpc[j + 1][i] = GLOB.tpc[2 * j + 1][i];
            }
        }
        *mdt *= 2;
        *idt = 6;
    }

    return 0;
}

// do .plg

int Temp::putemp(int iwn, double *tp, double tchcal, int iw, int nw, double tempgrd, double *dtempgrd)
{

    for (int i = 0; i < iwn; ++i)
    {
        GLOB.valxn[iwn + i] = GLOB.valxn[i];
        GLOB.valyn[iwn + i] = GLOB.valyn[i];
    }
    setval(tp, nw, iw, iwn);

    double dtmax = 0.;
    for (int i = 0; i < iwn; ++i)
    {
        dtmax = std::max(dtmax, std::abs(GLOB.valxn[i] - GLOB.valxn[iwn + i]));
    }
    *dtempgrd += dtmax;
    if (*dtempgrd >= tempgrd)
        //if ((*dtempgrd > tempgrd) || std::abs(*dtempgrd - tempgrd) < 1e-10)
    {
        *dtempgrd = 0.0;

        std::vector<double> v_temp;
        v_temp.push_back(tchcal);
        for (int j = 0; j < iwn; ++j)
        {
            v_temp.push_back(GLOB.valxn[j]);
        }
        v_plg.push_back(v_temp);

    }

    return 0;
}

// do .tlg

int Temp::putem(int istem, double *f, const std::vector<double>& tp, double tchcal, int *nstem, double temgrd, double *dtemgrd)
{

    int it;
    double dtmax = 0.;

    for (int i = 0; i < istem; ++i)
    {
        it = nstem[i];
        dtmax = std::max(dtmax, std::abs(f[it - 1] - tp[it - 1]));
    }

    *dtemgrd += dtmax;
    if (*dtemgrd >= temgrd)
        //if ((*dtemgrd > temgrd) || std::abs(*dtemgrd - temgrd) < 1e-10)
    {
        *dtemgrd = 0.0;

        std::vector<double> v_temp;
        v_temp.push_back(tchcal);
        for (int j = 0; j < istem; ++j)
        {
            v_temp.push_back(f[nstem[j] - 1]);
        }
        v_tlg.push_back(v_temp);
    }

    return 0;
}

int Temp::setval(double *tp, int nw, int iw, int iwn)
{

    int ix, iy;

    for (int i = 0; i < iw; ++i)
    {
        iy = 2 * (i + 1);
        ix = iy - 1;
        GLOB.x[i] = GLOB.zargd.xw[ix - 1];
        GLOB.y[i] = GLOB.zargd.xw[iy - 1];
        GLOB.valx[i] = tp[i];
        GLOB.valy[i] = 0.;
    }

    setvl(GLOB.x, GLOB.y, GLOB.xn, GLOB.yn, GLOB.valx,
          GLOB.valy, GLOB.valxn, GLOB.valyn, nw, iw, iwn);

    return 0;
}


/* ----> Podprogram przygotowujacy dane dla procedury okreslajacej */
/* ----> wartosc parametrow (valxn(i),valyn(i)) w i-tym wezle nowej sieci */
/* ----> elementow (xn(i),yn(i)) na podstawie ich wartosci (valx,valy) w */
/* ----> czterech najblizszych wezlach starej sieci (x,y). Ilosc wezlow */

/* ----> nowej sieci - iwn, a starej iw (w nw wierszach) */


int Temp::setvl(double *x, double *y, double *xn,
                double *yn, double *valx, double *valy, double *valxn,
                double *valyn, int nw, int iw, int iwn)
{

    double v[4];
    int ij;
    bool ok;
    double xx[4], yy[4], odl;
    int num[4];
    int numb[4];
    double dist;
    int inum, numb1;
    double dist1, odlmin, valmin, valmax;
    bool minmax = true;

    for (int i = 0; i < iwn; ++i)
    {
        numb[0] = 1;
        dist1 = 1e15;

        for (int j = 0; j < iw; ++j)
        {
            dist = std::sqrt(pow(xn[i] - x[j], 2) + pow(yn[i] - y[j], 2));
            if (dist < dist1)
            {
                dist1 = dist;
                numb[0] = j + 1;
            }
        }
        if (std::abs(dist1) < 1e-15)
        {
            valxn[i] = valx[numb[0] - 1];
            valyn[i] = valy[numb[0] - 1];
        }
        else
        {
            numb1 = numb[0];
            odlmin = 1e10;
            for (int j = 1; j <= 4; ++j)
            {
                if (j == 1)
                {
                    num[0] = numb1 - nw - 1;
                }
                if (j == 2)
                {
                    num[0] = numb1 - 1;
                }
                if (j == 3)
                {
                    num[0] = numb1;
                }
                if (j == 4)
                {
                    num[0] = numb1 - nw;
                }
                inum = num[0] / nw;
                inum = num[0] - inum * nw;
                ok = num[0] > 0 && num[0] < iw - nw && inum > 0;
                if (ok)
                {
                    num[1] = num[0] + nw;
                    num[2] = num[1] + 1;
                    num[3] = num[0] + 1;
                    odl = odlg(num, x, y, xn[i], yn[i]);
                    if (odl < odlmin)
                    {
                        odlmin = odl;
                        for (int k = 0; k < 4; ++k)
                        {
                            numb[k] = num[k];
                        }
                    }
                }
            }
            valmin = 1e10;
            valmax = -1e10;
            for (int j = 0; j < 4; ++j)
            {
                xx[j] = x[numb[j] - 1];
                yy[j] = y[numb[j] - 1];
                v[j] = valx[numb[j] - 1];
                if (v[j] < valmin)
                {
                    valmin = v[j];
                }
                if (v[j] > valmax)
                {
                    valmax = v[j];
                }
            }
            ij = 0;
            interp(xx, yy, v, &xn[i], &yn[i], &valxn[i], ij);
            if (minmax)
            {
                if (valxn[i] > valmax)
                {
                    valxn[i] = valmax;
                }
                if (valxn[i] < valmin)
                {
                    valxn[i] = valmin;
                }
            }
            valmin = 1e10;
            valmax = -1e10;
            for (int j = 0; j < 4; ++j)
            {
                v[j] = valy[numb[j] - 1];
                if (v[j] < valmin)
                {
                    valmin = v[j];
                }
                if (v[j] > valmax)
                {
                    valmax = v[j];
                }
            }
            ij = 1;
            interp(xx, yy, v, &xn[i], &yn[i], &valyn[i], ij);
            if (minmax)
            {
                if (valyn[i] > valmax)
                {
                    valyn[i] = valmax;
                }
                if (valyn[i] < valmin)
                {
                    valyn[i] = valmin;
                }
            }
        }
    }
    return 0;
}

double Temp::odlg(int *num, double *x, double *y, double x0, double y0)
{

    double dx, dy, odl = 0.0;

    for (int i = 0; i < 4; ++i)
    {
        dx = x[num[i] - 1] - x0;
        dy = y[num[i] - 1] - y0;
        odl += std::sqrt(dx * dx + dy * dy);
    }

    return odl;
}



/* ---->  Podprogram okresla wartosc uu dla danych wspolrzednych (xx,yy) */
/* ---->  na podstawie znanyxh wartosci u(i) w punktach (x(i),y(i)) */
/* ---->  otaczajacych punkt (xx,yy). Wskaznik ji = 0 powoduje wykonanie */
/* ---->  wstepnych obliczen w przypadku gdy uklad punktow otaczajacych */
/* ---->  sie zmienil. */
/* ---->   W przypadku gdy wyzn jest rowny 0. punkt okreslany jest na */

/* ---->   podstawie tylko dwu pierwszych punktow (x(i),y(i)). */

int Temp::interp(double *x, double *y, double *u, double *xx, double *yy, double *uu, int ij)
{

    if (ij == 0)
    {
        double x12, x13, x14, x23, x24, x34, y12, y13, y14, y23, y24,
                y34, z11, z12, z13, z21, z22, z23, z31, z32, z33, z41, z42, z43;

        x12 = x[0] - x[1];
        x13 = x[0] - x[2];
        x14 = x[0] - x[3];
        x23 = x[1] - x[2];
        x24 = x[1] - x[3];
        x34 = x[2] - x[3];
        y12 = y[0] * y[1] * x12;
        y13 = y[0] * y[2] * x13;
        y14 = y[0] * y[3] * x14;
        y23 = y[1] * y[2] * x23;
        y24 = y[1] * y[3] * x24;
        y34 = y[2] * y[3] * x34;
        static_.a[0] = x[1] * y34 - x[2] * y24 + x[3] * y23;
        static_.a[4] = -x[0] * y34 + x[2] * y14 - x[3] * y13;
        static_.a[8] = x[0] * y24 - x[1] * y14 + x[3] * y12;
        static_.a[12] = -x[0] * y23 + x[1] * y13 - x[2] * y12;
        static_.a[1] = -y23 + y24 - y34;
        static_.a[5] = y13 - y14 + y34;
        static_.a[9] = -y12 + y14 - y24;
        static_.a[13] = y12 - y13 + y23;
        z11 = y[0] * x23;
        z12 = y[0] * x24;
        z13 = y[0] * x34;
        z21 = y[1] * x13;
        z22 = y[1] * x14;
        z23 = y[1] * x34;
        z31 = y[2] * x12;
        z32 = y[2] * x14;
        z33 = y[2] * x24;
        z41 = y[3] * x12;
        z42 = y[3] * x13;
        z43 = y[3] * x23;
        static_.a[2] = x[1] * z23 - x[2] * z33 + x[3] * z43;
        static_.a[6] = -x[0] * z13 + x[2] * z32 - x[3] * z42;
        static_.a[10] = x[0] * z12 - x[1] * z22 + x[3] * z41;
        static_.a[14] = -x[0] * z11 + x[1] * z21 - x[2] * z31;
        static_.a[3] = -z23 + z33 - z43;
        static_.a[7] = z13 - z32 + z42;
        static_.a[11] = -z12 + z22 - z41;
        static_.a[15] = z11 - z21 + z31;
        static_.wyzn = static_.a[0] + static_.a[4] + static_.a[8] + static_.a[12];
    }
    //    if (static_.wyzn != 0.)
    if (std::abs(static_.wyzn) > epsilon)
    {
        double wsp[4];
        for (int j = 1; j <= 4; ++j)
        {
            wsp[j - 1] = 0.;
            for (int i = 1; i <= 4; ++i)
            {
                wsp[j - 1] += static_.a[j + 4 * i - 5] * u[i];
            }
            wsp[j - 1] /= static_.wyzn;
        }
        *uu = wsp[0] + wsp[1] * *xx + wsp[2] * *yy + wsp[3] * *xx * *yy;
    }
    else
    {
        double x1_2__ = x[0] - x[1];
        if (std::abs(x1_2__) < 1e-10)
        {
            *uu = (u[0] + u[1]) / 2.;
        }
        else
        {
            *uu = (u[0] * (*xx - x[1]) + u[1] * (x[0] - *xx)) / x1_2__;
        }
    }
    return 0;
}



/*   DT -  krok czasowy */
/*   CRATE - predkosc chlodzenia */
/*   TIME - czas biezacy */
/*   TTRA - temperatura przemiany */
/*   TIP0 - czas poczatku przemiany (TABLICA) */
/*   XV   - ulamek przemieniony */
/*   Q    - cieplo przemiany w kroku dt */
/*   DZ   - pocz. wielkosc ziarna */
/*   Tau  - czas inkubacji przemiany perlitycznej */

/*          (dla modelu addytywnego) */

double Temp::cieprzem(double dt, double t, double crate, double time, double *tip0, double *xvv)
{

    double dz = 50.;
    double b, q, x, y, dh, xv1, t1000, rho, ttra;

    GLOB.przemxv = *xvv;
    if (*tip0 < 0.)
    {
        // nieuzywane !!!!!!!!1
        //double xmn = 1.2;
        //double xk = xmn * -.8 * xmn + xmn * 1.86 - 1.04;
        //double t0 = xmn * -225. * xmn + xmn * 452.5 + 746.;
        //double tau = xk * t * exp(9812. / t) / pow(t0 - t, 2.8);

        /*       dla 1.1% Mn */

        /*            TtRa=645.8*CRate**(-0.0165) */

        /*       dla 1.2% Mn */

        ttra = pow(crate, -.0165) * 639.6;


        /*            TTRA=631.9*CRATE**(-.02434) */
        if (ttra > 800.)
        {
            ttra = 800.;
        }
        if (ttra > t)
        {
            *tip0 = time - dt;
        }
    }

    if (*tip0 < 0.)
    {
        /*         IF(Tip0.lt.0.d0 .or. Time.LT.Tip0+Tau)  then */
        q = 0.;
    }
    else
    {
        /*            b=-3.26d-4*T*T+.361*T-95.37 */
        b = t * -2.55e-4 * t + t * .2794 - 72.23;
        b = exp(b);
        b *= .5;
        //d__1 = time - *tip0;
        xv1 = 1. - exp(-b * pow(time - *tip0, 2.04) / pow(dz, 1.96));
        if (xv1 < GLOB.przemxv)
        {
            xv1 = GLOB.przemxv;
        }
        dh = 245.16 - (t + 273.) * .168;
        t1000 = t / 1000.;
        alrola(t1000, &x, &rho, &y);
        q = rho * dh * (xv1 - GLOB.przemxv) / dt;
        GLOB.przemxv = xv1;
        q *= 1000.;
    }
    *xvv = GLOB.przemxv;

    return q;
}


// iloczyn
// a(n,nm) * b(nm,m) = c(n,m)

int Temp::ilo(double *a, double *b, int n, int nm, int m, double *c)
{
    /* float *a = new float[nm * n];
     for (int i = 0; i < n*nm; ++i)
         a[i] = static_cast<float> (aa[i]);
         //for (int k = 0; k < nm; ++k)
           //  a[i + k * n] = static_cast<float> (aa[i + k * n]);

     float *b = new float[nm * m];
     for (int j = 0; j < nm*m; ++j)
        // for (int k = 0; k < nm; ++k)
             //b[j + k * n] = static_cast<float> (bb[j + k * n]);
         b[j] = static_cast<float> (bb[j]);
     */
    //float *c = new float[n+m*n];

    /* cl::Buffer A(*GPU.context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof (cl_float) * (nm * n), a);
     cl::Buffer B(*GPU.context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof (cl_float) * (nm * m), b);
     cl::Buffer C(*GPU.context, CL_MEM_WRITE_ONLY | CL_MEM_USE_HOST_PTR, sizeof (cl_float) * (m * n), c);
     //GPU.queue->finish();

     cl::Kernel k(*GPU.program, "ilo");
     k.setArg(0, A);
     k.setArg(1, B);
     k.setArg(2, n);
     k.setArg(3, nm);
     k.setArg(4, m);
     k.setArg(5, C);

   //  GPU.queue->enqueueWriteBuffer(A, CL_TRUE, 0, sizeof (cl_float) * (n + nm * n), a);
   //  GPU.queue->enqueueWriteBuffer(B, CL_TRUE, 0, sizeof (cl_float) * (nm + nm * m), b);

     // GPU.queue->finish();
     //const size_t global_work_size[2] = {n + m*n,n + m*n};
    //const size_t local_work_size[2] = {16,16};
    //GPU.queue->enqueueNDRangeKernel(k, cl::NullRange, cl::NDRange(n + m*n), cl::NDRange(2, 2));
     GPU.queue->enqueueNDRangeKernel(k, cl::NullRange, cl::NDRange(4,4), cl::NDRange(4,4));
    //GPU.queue->enqueueNDRangeKernel(k, cl::NullRange, cl::NDRange(nm), cl::NDRange(2,2));


     GPU.queue->enqueueReadBuffer(C, true, 0, sizeof (cl_float) * (m * n), c);
     //  GPU.queue->finish();
     //char *p = (char*)GPU.queue->enqueueMapBuffer(buffer1, true, CL_MAP_READ, 0, 16);
     //std::cout << "p = " << p << std::endl;
     //GPU.queue->enqueueUnmapMemObject(buffer1, p);
      //
     */


    for (int i = 0; i < n; ++i)
    {
        for (int j = 0; j < m; ++j)
        {
            c[i + j * n] = 0.;
            for (int k = 0; k < nm; ++k)
            {
                c[i + j * n] += a[i + k * n] * b[k + j * nm];
            }
        }
    }

    // std::cout << "OUT " << static_cast<float>(c[0]) << std::endl;

    // delete[] a;
    // delete[] b;
    return 0;
}

int Temp::elem(bool iosym, int nel, int it, int *iw, double rlchar, double *dysyp, double dt, int *npas1, int *ij, int *ji, int qq)
{

    int iciep = 0;
    int iwws = 0;
    double ksi[4] = {-.57735026919, .57735026919, .57735026919, -.57735026919};
    double eta[4] = {-.57735026919, -.57735026919, .57735026919, .57735026919};
    double xksi[2] = {-.57735026919, .57735026919};
    //double ksi[4] = {-.57735, .57735, .57735, -.57735};
    //double eta[4] = {-.57735, -.57735, .57735, .57735};
    //double xksi[2] = {-.57735, .57735};
    int i, j, l, m;
    double n[4], s[4] /* bylo [2][2] */, x;
    int i1, j1;
    double aa, ca, bb, cg[16] /* bylo [4][4] */, ei[1] /* 
	    bylo [1][1] */, cm[16] /* bylo [4][4] */, dl;
    int ii[2];
    double hm[16] /* bylo [4][4] */, tb[2], fm[4];
    int il;
    double bx[4], by[4], nn[2], dx, dy, sm[16] /* bylo [4][4] 
	    */, hx[16] /* bylo [4][4] */, hy[16] /* bylo [4][4] */, ss,
            tt[4], qw[4], rx, qp, xx[4], sa, ro, cp, dt3;
    int ibb[4];
    double eim, tei[1] /* bylo [1][1] */, vei[1] /* 
	    bylo [1][1] */;

    double ste[1] /* bylo [1][1] */, ssi;
    double tsr[1] /* bylo [1][1] */, xww[8];
    double detj;
    int kkod[4];
    double ndet[4], teim, epsi[4], veim, stem, xvel[4], tsrm;
    double vepsi[4], razpr;
    double lambda;

    dt3 = 3. / dt;
    for (i = 1; i <= 4; ++i)
    {
        xww[i * 2 - 2] = GLOB.zargd.xw[GLOB.iy[nel - 1][i - 1] * 2 - 2];
        xww[i * 2 - 1] = GLOB.zargd.xw[GLOB.iy[nel - 1][i - 1] * 2 - 1];
        ibb[i - 1] = GLOB.ib[GLOB.iy[nel - 1][i - 1] - 1];
        tt[i - 1] = GLOB.tp[GLOB.iy[nel - 1][i - 1] - 1];
        epsi[i - 1] = GLOB.tweps[GLOB.iy[nel - 1][i - 1] - 1];
        vepsi[i - 1] = GLOB.dteps[GLOB.iy[nel - 1][i - 1] - 1];
        qw[i - 1] = GLOB.q[GLOB.iy[nel - 1][i - 1] - 1];
        xvel[i - 1] = GLOB.xvv[GLOB.iy[nel - 1][i - 1] - 1];
    }
    /*     if(it.eq.1.or.ralfpo.gt.5.0) then */
    if (it == 1 || GLOB.ral.ralfpo[qq] > .5)
    {
        ss = 0.;
        for (i = 1; i <= 4; ++i)
        {
            fm[i - 1] = 0.;
            ndet[i - 1] = 0.;
            for (j = 1; j <= 4; ++j)
            {
                hm[i + 4 * j - 5] = 0.;
                cm[i + 4 * j - 5] = 0.;
                sm[i + 4 * j - 5] = 0.;
            }
        }
        for (l = 1; l <= 4; ++l)
        {
            fmatb(xww, ksi[l - 1], eta[l - 1], bx, by, n, &detj);

            /* ---->  stan osiowo-symetryczny */

            if (iosym)
            {
                rx = 0.;
                for (i = 1; i <= 4; ++i)
                {
                    rx += n[i - 1] * xww[2 * i - 2];
                }
                detj = rx * detj;
            }
            ilo(n, epsi, 1, 4, 1, ei);
            eim = ei[0];
            ilo(n, vepsi, 1, 4, 1, vei);
            veim = vei[0];
            ilo(n, tt, 1, 4, 1, tei);
            teim = tei[0];
            ilo(n, qw, 1, 4, 1, tei);
            qp = tei[0];
            ilo(n, xvel, 1, 4, 1, tei);
            GLOB.przemxv = tei[0];
            /*           if (ralfpo.lt.5.d0) then */
            ssi = 0.;
            /*           else */
            /*              ssi= sigma(eim,veim,teim)*snorm */
            /*              ssi= dysyp*1.d6*ssi */
            /*           end if */
            ilo(bx, bx, 4, 1, 4, hx);
            ilo(by, by, 4, 1, 4, hy);
            ilo(n, n, 4, 1, 4, cg);
            ss += detj;
            for (i = 1; i <= 4; ++i)
            {
                ndet[i - 1] += n[i - 1] * detj;
                /*              fm(i)=fm(i)-ssi*veim*n(i)*detj */
                fm[i - 1] -= qp * n[i - 1] * detj;
                for (j = 1; j <= 4; ++j)
                {
                    hm[i + 4 * j - 5] += (hx[i + 4 * j - 5] + hy[i + 4 * j - 5]) * detj;
                    cm[i + 4 * j - 5] += cg[i + 4 * j - 5] * detj;
                }
            }
        }
        ilo(ndet, tt, 1, 4, 1, ste);
        stem = ste[0];
        ++(*ij);
        cnd[*ij - 1] = ss;
        for (i = 0; i < 4; ++i)
        {
            ++(*ij);
            cnd[*ij - 1] = ndet[i];
            /* L401: */
        }
        for (i = 1; i <= 4; ++i)
        {
            for (j = i; j <= 4; ++j)
            {
                *ji += 2;
                hcf[*ji - 2] = hm[i + 4 * j - 5];
                hcf[*ji - 1] = cm[i + 4 * j - 5];
            }
            ++(*ji);
            hcf[*ji - 1] = fm[i - 1];
        }
    }
    else
    {
        for (i = 1; i <= 4; ++i)
        {
            for (j = 1; j <= 4; ++j)
            {
                sm[i + (j << 2) - 5] = 0.;
            }
        }
        for (i = 1; i <= 4; ++i)
        {
            for (j = i; j <= 4; ++j)
            {
                *ji += 2;
                hm[i + (j << 2) - 5] = hcf[*ji - 2];
                cm[i + (j << 2) - 5] = hcf[*ji - 1];
                hm[j + (i << 2) - 5] = hm[i + (j << 2) - 5];
                cm[j + (i << 2) - 5] = cm[i + (j << 2) - 5];
            }
            ++(*ji);
            fm[i - 1] = hcf[*ji - 1];
        }
        ++(*ij);
        ss = cnd[*ij - 1];
        for (i = 1; i <= 4; ++i)
        {
            ++(*ij);
            ndet[i - 1] = cnd[*ij - 1];
        }
        ilo(ndet, tt, 1, 4, 1, ste);
        stem = ste[0];
    }
    if (GLOB.kod[nel - 1] != 0)
    {
        dkod(GLOB.kod[nel - 1], &il, kkod);
        if (iwws == 1 && it == 1)
        {
            FM.FL2 << nel << "\n";
        }

        for (i = 1; i <= il; ++i)
        {
            ii[0] = kkod[i - 1];
            ii[1] = ii[0] + 1;
            if (ii[1] == 5)
            {
                ii[1] = 1;
            }
            xx[0] = xww[(ii[0] << 1) - 2];
            xx[1] = xww[(ii[0] << 1) - 1];
            xx[2] = xww[(ii[1] << 1) - 2];
            xx[3] = xww[(ii[1] << 1) - 1];
            tb[0] = tt[ii[0] - 1];
            tb[1] = tt[ii[1] - 1];
            dx = xx[2] - xx[0];
            dy = xx[3] - xx[1];
            dl = std::sqrt(dx * dx + dy * dy);

            /* ---->  Sprawdzenie czy brzeg elementu styka sie z narzedziem */

            if (ibb[ii[0] - 1] >= 10 && ibb[ii[1] - 1] >= 10)
            {
                /*              if(dx.lt.dy) then */
                /*                 razpr=ralfpn */
                /*              else */
                /*                 razpr=ralfpo */
                /*              end if */
                /*           else if (ralfa.gt.0.d0) then */
                /*              razpr=1.d0 */
                razpr = GLOB.ral.ralfa[qq];
            }
            else
            {
                razpr = 0.;
            }
            GLOB.ral.ralfae = razpr;
            if (iwws == 1 && it == 1)
            {
                if (GLOB.ral.ralfae > 1.)
                {
                    FM.FL2 << "chlodzenie przez walec na odcinku miedzy wezlami " << ii[0] << " i " << ii[1] << "\n";
                }
                else
                {
                    FM.FL2 << "chlodzenie powietrzem na odcinku miedzy wezlami " << ii[0] << " i " << ii[1] << "\n";
                }
            }
            sa = dy / dl;
            ca = dx / dl;
            rx = (xx[2] + xx[0]) / 2.;
            GLOB.rot(xx, &xx[1], sa, ca);
            GLOB.rot(&xx[2], &xx[3], sa, ca);
            dx = std::abs(xx[2] - xx[0]);
            /*           detj=(xx(3)-xx(1))/2. */
            detj = std::abs((xx[2] - xx[0]) / 2.);

            /* ---->  stan osiowo-symetryczny */

            if (iosym)
            {
                detj = rx * detj;
            }
            aa = (xx[2] - xx[0]) / 2.;
            bb = (xx[2] + xx[0]) / 2.;
            for (j = 0; j < 2; ++j)
            {
                x = aa * xksi[j] + bb;
                nn[0] = (xx[2] - x) / dx;
                nn[1] = (x - xx[0]) / dx;
                ilo(nn, nn, 2, 1, 2, s);
                ilo(nn, tb, 1, 2, 1, tsr);
                tsrm = tsr[0];

                double alfazp;
                double alfaz;
                if (GLOB.ral.ispray)
                {
                    alspray(tsrm, &rlchar, &razpr, &alfaz, &alfazp);
                }
                else
                {
                    alrol1(&tsrm, rlchar, razpr, &alfaz, &alfazp);
                }
                if (iciep == 1)
                {
                    FM.FL2 << "element nr. " << nel << " tsr= " << tsrm << "oC ,alfa= " << alfaz << " J/Km^2\n";
                }
                for (l = 1; l <= 2; ++l)
                {
                    fm[ii[l - 1] - 1] -= nn[l - 1] * detj * alfazp;
                    for (m = 1; m <= 2; ++m)
                    {
                        sm[ii[l - 1] + (ii[m - 1] << 2) - 5] += s[l + (m << 1) - 3] * detj * alfaz;
                    }
                }
            }
        }
    }
    tsrm = stem / ss;

    alrola(tsrm / 1000., &lambda, &ro, &cp);

    if (iciep == 1)
    {
        FM.FL2 << "element nr. " << nel << " tsr= " << tsrm << "oC     lambda= " << lambda
                << " W/mK    ro= " << ro << " kg/m^3  cp=  " << cp << " J/kgK  \n";
    }
    for (i = 1; i <= 4; ++i)
    {
        for (j = 1; j <= 4; ++j)
        {
            hm[i + (j << 2) - 5] = hm[i + (j << 2) - 5] * lambda + sm[i + (j << 2) - 5];
            cm[i + (j << 2) - 5] *= ro * cp;
        }
    }
    for (i = 1; i <= 4; ++i)
    {
        i1 = GLOB.iy[nel - 1][i - 1];
        f[i1 - 1] += fm[i - 1];
        for (j = 1; j <= 4; ++j)
        {
            j1 = GLOB.iy[nel - 1][j - 1];
            if (j1 >= i1)
            {
                j1 += -i1 + 1;
                *npas1 = std::max<int>(*npas1, j1);
                //                c___.c[j1 - 1][i1 - 1] += dt3 * cm[i + (j << 2) - 5] - hm[i + (j << 2) - 5];
                //                h___.h[j1 - 1][i1 - 1] += dt3 * cm[i + (j << 2) - 5] + hm[i + (j << 2) - 5] * 2.;
                c[i1 + j1 * MWMAX * NWMAX - (MWMAX * NWMAX + 1)] += dt3 * cm[i + (j << 2) - 5] - hm[i + (j << 2) - 5];
                h[i1 + j1 * MWMAX * NWMAX - (MWMAX * NWMAX + 1)] += dt3 * cm[i + (j << 2) - 5] + hm[i + (j << 2) - 5] * 2.;
            }
        }
    }
    return 0;
}

int Temp::dkod(int kod, int *ikod, int *kkod)
{
    int k0 = kod, k1;

    *ikod = 0;

    while (true)
    {
        k1 = k0 / 10;
        ++(*ikod);
        kkod[*ikod - 1] = k0 - k1 * 10;
        if (kkod[*ikod - 1] == 0)
        {
            --(*ikod);
        }
        k0 = k1;
        if (k0 == 0)
        {
            return 0;
        }
    }
}

int Temp::fn(double ksi, double eta, double *n)
{
    double ksii[4] = {-1., 1., 1., -1.};
    double etai[4] = {-1., -1., 1., 1.};

    for (int i = 0; i < 4; ++i)
    {
        n[i] = (ksii[i] * ksi + 1.) * (etai[i] * eta + 1.) / 4.;
    }

    return 0;
}



/* ----> podprogram wyznacza dla danego elementu i dla danego punktu */
/* ----> (ksi,eta) w lokalnym ukladzie wspolrzednych wartosci jakobianu */
/* ----> [B] a takze wartosci funkcji ksztaltu {N} dla temperatur. */

/* ----> xww - wektor wspolrzednych koncowych (parametrow) dla elementu */

int Temp::fmatb(double *xww, double ksi, double eta,
                double *bx, double *by, double *n, double *detj)
{

    double x12 = xww[0] - xww[2];
    double x13 = xww[0] - xww[4];
    double x41 = xww[6] - xww[0];
    double x23 = xww[2] - xww[4];
    double x24 = xww[2] - xww[6];
    double x34 = xww[4] - xww[6];
    double y12 = xww[1] - xww[3];
    double y13 = xww[1] - xww[5];
    double y41 = xww[7] - xww[1];
    double y23 = xww[3] - xww[5];
    double y24 = xww[3] - xww[7];
    double y34 = xww[5] - xww[7];
    double vol = x13 * y24 - x24 * y13;

    *detj = std::abs(vol + (x34 * y12 - x12 * y34) * ksi + (x41 * y23 - x23 * y41) * eta);
    bx[0] = (y24 - y34 * ksi - y23 * eta) / *detj;
    bx[1] = (y34 * ksi - y13 - y41 * eta) / *detj;
    bx[2] = (y12 * ksi - y24 + y41 * eta) / *detj;
    bx[3] = (y13 - y12 * ksi + y23 * eta) / *detj;
    by[0] = (x34 * ksi - x24 + x23 * eta) / *detj;
    by[1] = (x13 - x34 * ksi + x41 * eta) / *detj;
    by[2] = (x24 - x12 * ksi - x41 * eta) / *detj;
    by[3] = (x12 * ksi - x13 - x23 * eta) / *detj;
    *detj /= 8.;

    fn(ksi, eta, n);

    return 0;
}

void Temp::AddT(const std::vector<double>& t)
{
    v_t.push_back(t);
}

void Temp::AddTFinal()
{
    v_tf.push_back(v_t.back());
}

void Temp::AddTchcal(double val)
{
    v_tchcal.push_back(val);
}

void Temp::ShowTemps(bool val)
{
    v_show_temps.push_back(val);
}


/*     ----------------------------------------- */
/*     TSR - temperatura/1000 */
/*     LAMBDA - wspolczynnik przewodzenia ciepla */
/*     RO - gestosc */
/*     CP -cieplo wlasciwe */
/*     ----------------------------------------- */

/*    1.084D0,-24.6D0,23.16D0,51.96D0,-2.02519D0/ */

/*    g - austenit, a - perlit */

int Temp::alrola(double tsr, double *lambda, double *ro, double *cp)
{
    /* Initialized data */

    double ros = 7850.;
    double cc1 = 207.9;
    double ca1 = 249.4;
    double cb1 = 1.1412;
    double cl = 17.222;
    double al = 61.923;
    double bl = -4.543;

    double tk = tsr * 1000. + 273.;
    /*         CPG=CC2+CA2*TSR**CB2 */
    double cpg = exp(tsr * .00378) * 46.16 + 626.2;
    double cpa = cc1 + ca1 * exp(cb1 * tsr);
    *cp = cpg * (1. - GLOB.przemxv) + cpa * GLOB.przemxv;
    /* -------- Mirek */
    double xla = cl + al * exp(bl * tsr);
    /* -------- Morales,Lopez,Olivares, ISIJ 1990,p.48 */
    double xlg = tk * .013 + 9.84;
    *lambda = xlg * (1. - GLOB.przemxv) + xla * GLOB.przemxv;
    /* Computing 3rd power */
    double d__1 = tsr * .008 * tsr / 2. + 1.;
    *ro = ros / (d__1 * (d__1 * d__1));

    /*     Nowe dane Romana ( osobno dla austenitu i perlitu ) */

    /*     RoG=8024.1-0.524*Tk */
    /*     RoA=7863.1-0.3118*Tk-6.0894d-5*Tk*Tk */
    /*     Ro=RoG*(1.-XV)+RoA*XV */
    return 0;
}


/*        alfa dla chlodzenia natryskiem */
/*        zwspr - strumien wody l/s/m/m */
/*  t     temperatura w stopniach Celsjusz */

/*  alfa  wspolczynnik wymiany ciepla W/K/m/m */


int Temp::alspray(double t, double *l, double *razpr, double *alfa, double *alfa1)
{

    double zwspr = 4.;
    double d__1 = (t - 700.) / (exp(t * .1 - 10.) + 1.) + 700.;
    double zmul = 1. / pow(d__1, 2.455);
    *alfa = pow(zwspr, .616) * 3.15e6 * zmul * (1. - 1. / (exp(t * .025 - 6.25) + 1.));

    /*        1.3 zwiekszenie efektywnosci chlodzenia przez dodanie powietrza */

    *alfa = *alfa * 1000. * 1.3;
    *alfa1 = *alfa * GLOB.ral.tpow;
    return 0;
}

int Temp::alrol1(double *tsr, double l, double razpr, double *alfa, double *alfa1)
{

    if (razpr > 5.)
    {
        *alfa = razpr;
    }
    else
    {
        double cs = 122.;
        double eta0 = 1.719e-5;
        double rop = 1.293;
        double beta = .003661;
        double g = 9.81;
        double eps = .8;
        double lmbdap = .0255;
        double cpp = 1005.;
        double d__1, d__2;

        double dtem = *tsr - GLOB.ral.tpow;
        if (dtem < .001)
        {
            dtem = .001;
            *tsr = GLOB.ral.tpow + dtem;
        }
        double t = *tsr + 273.;
        double eta = eta0 * std::sqrt(t / 273.) * (cs / 273. + 1.) / (cs / t + 1.);
        double gr = g * l * l * l * rop * rop * beta * dtem / eta / eta;
        double pr = cpp * eta / lmbdap;
        double grpr = gr * pr;

        if (grpr <= 500.)
        {
            double d__3;
            d__2 = l, d__3 = d__2, d__2 *= d__2;
            d__1 = dtem / (d__3 * (d__2 * d__2));
            double alf = exp(*tsr * 3.416e-4) * .3022 * pow(d__1, .125);
        }
        if (grpr > 500. && grpr <= 2e7)
        {
            double alf = (exp(*tsr * -7.74e-4) * 1.0804 + .3211) * pow(dtem / l, .25);
        }
        if (grpr > 2e7)
        {
            double alf = (exp(*tsr * -.002844) * 1.5745 + .4958) * pow(dtem, .33333);
        }

        /*       wspolczynniki emisji wg. Zbyszka */
        eps = *tsr / 1000. * (*tsr * .125 / 1000. - .38) + 1.1;

        /* Computing 4th power */
        d__1 = t;
        d__1 *= d__1;
        /* Computing 4th power */
        d__2 = GLOB.ral.tpow + 273., d__2 *= d__2;
        double alfpr = eps * 5.675e-8 * (d__1 * d__1 - d__2 * d__2) / dtem;
        alfpr *= razpr;
        /*       alfa=(alf+alfpr)*ralfae */
        *alfa = alfpr * GLOB.ral.ralfae;
    }
    *alfa1 = *alfa * GLOB.ral.tpow;
    return 0;
}
